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SECTION  1 


INTRODUCTION 

In  an  effort  to  aid  in  the  design  and  fabrication  of  devices  more  resistant  to 
single  event  upsets  and  to  gain  understanding  of  the  internal  dynamics  of 
devices  struck  by  single  radiation  particles,  device  researchers  have  turned 
to  numerical  simulation.  Early  studies  of  such  phenomena  involved 
two-dimensional  simulations  of  the  response  of  two-terminal  N+P  diode 
structures  to  single  particle  radiation  -{1,2,3],  These  studies  gave  light  to 
a  result  coined  the  field- funneling  effect.  ' 

However,  the  picture  which  emerged  from  the  initial  diode  simulations  was 
incomplete.  As  was  demonstrated  in  [A, 5]  the  response  of  more  complex  devices 
was  significantly  different  from  that  of  simple  diodes.  For  example,  in  2-D 
simulations  of  JFETS  and  MESFETS  [4]  it  was  found  that  the  spreading  of  excess 
charge  deposited  in  the  device  substrate  lowered  the  substrate  resistance 
substantially  and  provided  a  source -drain  current  path  in  a  device  that  was 
initially  in  the  off  state.  Additionally,  only  a  small  portion  of  the  charge 
deposited  by  the  radiation  particle  was  collected  at  a  struck  gate  node. 
Similarly,  simulations  of  CMOS  and  bipolar  or  multi-layered  devices  [4,5] 
showed  the  existence  of  extremely  complex  current  paths,  aside  from  those  at 
the  struck  node,  and  including  such  phenomena  as  plama  wires  or  ion  shunts 
across  the  base  region  of  bipolar  devices.  Such  effects  have,  as  a  result, 
been  included  in  circuit  simulations  [6], 

Even  with  the  additional  information  provided  by  2-D  device  simulation, 
modeling  of  SEUs  using  circuit  simulations  is  limited.  There  is  substantial 
interaction  between  the  external  circuit  and  the  device.  Thus,  it  is 
important  to  couple  the  device  simulation  to  the  external  circuit.  Such  an 
approach  has  been  followed  in  [7]  where  a  two-dimensional  simulation  of  a  CMOS 
memory  cell,  coupled  to  the  external  circuit,  was  performed. 

While  the  present  authors  agree  strongly  that  drift  and  diffusion  simulations 
of  the  devices  in  question  must  be  coupled  to  the  external  circuitry,  the 
approach  taken  in  [7]  still  has  a  major  shortcoming;  the  drift  and  diffusion 
simulations  were  limited  to  two  dimensions.  While  it  is  true  that  much  has 
been  learned  through  two-dimensional  simulations  about  device  response 
characteristics  when  subjected  to  single  particle  radiation,  the  physics 
involved  is  inherently  three-dimensional.  In  [8,9]  the  present  authors 
demonstrated  the  significant  inaccuracies  resulting  from  the  assumption  of 
two-dimensionality  by  carefully  comparing  results  from  both  two-  and 
three-dimensional  simulations  of  charge  collection  in  a  simple  silicon  diode. 
It  was  also  demonstrated  that  the  two-dimensional  results  for  the  collection 
time  and  the  current  pulse  could  be  made  to  follow  the  three-dimensional 
result  by  properly  scaling  the  track  density.  Such  a  scaling  result  was 
recently  found  to  apply  to  an  NMOS  structure  which  was  simulated  in  both  two 
and  three  dimensions  in  [10].  However,  even  in  these  cases  it  was  cautioned 
that  the  transport  process  in  two  dimensions  was  highly  inaccurate.  Thus, 
while  two-dimensional  simulations  may  provide  a  useful  tool,  results  thus 
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obtained  must  be,  at  least,  verified  through  selective  three-dimensional 
computations . 

While  the  need  to  perform  three-dimensional  transient  simulations  is  readily 
apparent,  a  major  drawback  has  been  the  computational  intensity  of  such 
simulations.  In  1981,  Buturla,  et.  al.  [11]  reported  on  a  three-dimensional 
transient  simulation  involving  the  reverse  recovery  of  an  ellipsoidal 
junction.  The  computational  segment  of  the  device  simulated  was  1.25  x  1.25  x 
2.5/im  and  a  finite  element  grid  of  600  nodes  with  810  elements  was  used. 

Thirty  time  steps  were  taken  and  the  simulation  required  approximately  5  hours 
of  CPU  time  [11]  on  an  IBM  370/168  computer.  By  1985,  the  situation  had 
improved.  The  three-dimensional  transient  simulation  of  charge  collection 
performed  by  the  present  authors  [8,9]  only  required  approximately  3.5  hours 
of  CPU  time  on  a  Cray  1  computer  using  a  finite  difference  procedure  with 
17,500  grid  points  and  350  time  steps.  The  algorithm  used  in  those 
simulations,  while  highly  efficient,  did  not  take  advantage  of  the  Cray's 
vector  architecture.  However,  the  algorithm  was  ideally  suited  for 
vectorization  and  it  was  estimated  that  when  vectorized,  run  times  would  be 
reduced  by  at  least  an  order  of  magnitude  making  practical  three-dimensional 
transient  simulations  a  reality  for  more  complex  structures. 

This  was  indeed  found  to  be  a  conservative  estimate  in  [10]  where  a  vectorized 
version  of  this  algorithm  was  applied  to  the  simulation  of  a  silicon  NMOS 
device.  The  simulation  of  the  NMOS  structure  utilized  a  three-dimensional 
mesh  of  21560  points.  1050  times  steps  were  taken  during  the  simulation  which 
required  only  1.44  hours  of  CPU  time  on  a  Cray  1  computer.  This  was  roughly 
16  times  faster  than  the  run  time  projected  using  a  scalar  version  of  the  same 
solution  algorithm. 

The  present  report  continues  to  address  the  need  for  three-dimensional 
simulations  of  single  particle  effects  and  two-dimensional  scaling  approaches 
through  the  application  of  this  vector  algorithm  to  the  simulation  of  a  GaAs 
J  FET  s  true  tur e . 
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SECTION  2 


ANALYSIS 


2.1  GOVERNING  EQUATIONS. 


The  discussion  which  follows  is  general  in  that  allowance  is  made  for  hetero¬ 
junctions.  However,  in  the  JFET  structure  considered,  no  he teroj unctions  are 
present.  The  governing  continuity  equations  thus  take  the  form 


3N 

at 


+  G  -  R 


(1) 


3P 
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+  G  -  R 


(2) 


where  the  current  densities  are  given  by 
Jn  =  effyin  V (V’+Fn)  +  DnVNj 

Jp  =  -e[ppp  V(V»+Fp)  +  DpVpj 


(3) 

(4) 


The  quantities  Fn  and  Fp  are  introduced  to  account  for  variations  in  the 
conduction  and  valence  oand  energy  levels,  and  are  related  to  the  electron 
affinity,  the  density  of  states  and  the  energy  band  gap  as  [12] 


Fn  “  h  (X  +  kT  lnNc) 


(5) 


v 


FP  = 


(X  +  Eg  -  kT  lnNv) 


(6) 


The  gradients  of  Fn  and  Fp  give  rise  to  local  "effective  fields"  at 
material  interfaces  which  may  augment  or  retard  drift  transport  across  the 
interface.  G  and  R  in  Eqs .  (1)  and  (2)  are  generation  and  recombination 
terms.  The  recombination  term  is  given  as 


R  = 


NP  -  Ni2 

r  p ( N+N  t )  +  rn(P+Ni) 


+  r (N+P) (NP-Ni2 ) 


(7) 


Here  rn  and  rp  are  recombination  lifetimes  whereas  r  is  a  recombination 


rate  constant.  The  first  and  second  terms  in  Eq.  (7)  represent 
Shockley-Read-Hall  and  Auger  recombination,  respectively.  For  GaAs,  the 
carrier  lifetimes  are  on  the  order  of  1  nsec. 


Generation  of  carriers  results  from  two  effects  In  the  present  simulations; 
impact  ionization,  and  generation  due  to  the  energy  absorbed  from  the  incident 
radiation  particle.  The  authors  have  shown  [8]  that  at  bias  levels  of  a  few 
volts  impact  ionization  in  a  simple  diode  would  not  contribute  to  excess 
carrier  generation.  Since  the  bias  levels  imposed  here  on  reversed  junctions 
are,  at  a  maximum,  on  the  order  of  one  volt,  reverse  breakdown  of  such 
junctions  is  not  anticipated  and  impact  ionization  should  not  effect  the 
results.  However,  this  process  is  modeled  in  the  present  simulations 
following  the  standard  approach  for  completeness.  Generation  due  to  impact 
ionization  is  expressed  as 

Gi  -  \  («nM'nl  +  «p|j'p|)  (8) 


where  a^,  (k  =  n,p)  is  given  by 


(9) 


The  values  of  the  constants  A,  b  and  m  for  GaAs  are  given  in  Table  1.  It  is 
noted  that  the  current  densities  used  in  Eq.  (8)  are  limited  to  the  drift 
component  of  the  particle  currents  only,  following  Sze  [13],  rather  than  the 
full  particle  current  (drift  plus  diffusion) .  The  present  authors  believe 
that  this  approach  is  more  realistic  because,  for  example,  in  a  reversed  bias 
diode  ideally  the  current  flow  is  zero  with  drift  and  diffusion  currents 
balancing  exactly.  As  a  result,  breakdown  of  the  diode  at  high  reverse  bias 
would  not  be  predicted.  When  only  the  drift  component  is  considered, 
avalanche  generation  will  occur  at  high  fields  and  reverse  breakdown  will 
occur . 


Generation  due  to  thermalization  by  the  incident  particle  is  modeled  in  a 
straight  forward  manner  as 
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Na(r) 

— -  exp  [— t(r)/ra) 


(10) 


Here  Na(r)  is  the  concentration  of  particles  generated  within  the  track, 
and  can  be  a  function  of  distance  along  the  track.  ra  represents  a  time 
constant  thermalization,  typically  taken  as  3  psec ,  and  t(r)  represents  the 
time  elapsed  from  when  the  radiation  particle  penetrates  the  device  to  a 
specific  point  along  its  track.  For  example,  a  5Mev  alpha  particle  possesses 
a  range  of  18.5jim  in  GaAs  and  travels  at  an  average  velocity  of 
approximately  1.56  x  109  cm/sec.  Thus,  on  average,  carrier  generation  will 
begin  at  the  entry  point  of  the  particle  track  1.18psec  before  generation  at 
the  end  of  the  track. 
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While  the  model  is  sufficiently  general  to  include  the  above-mentioned 
effects,  in  the  present  simulations  NQ  was  assumed  constant  along  the 
entire  track  length,  and  the  delay  time  along  the  track  was  neglected. 
Therefore,  generation  was  assumed  to  take  place  simultaneously  and  uniformly 
along  the  entire  track  length. 

Since  space  charge  effects  must  be  considered  in  the  present  analysis,  a 
self-consistent  electric  field  must  be  determined  from  Poisson's  equation 

V*eV0  -  -p  -  e(N-P-C)  (U) 

where  C  is  the  net  doping  distribution  of  donor  and  acceptor  ions.  Here  we 
note  that,  due  to  the  hetero junction  formulation,  the  permittivity  could  be 
spatially  dependent. 


2.2  MOBILITY  AND  DIFFUSIVITY  MODELS. 


The  mobility  model  used  in  the  present  simulations  allows  for  negative 
differential  mobility  of  GaAs  electrons.  The  electron  velocity  is  given  by 


/n  -  [MoE  +  aCE/E^  +  ME/E^3  +  c(E/Ev) 4  J  —  (E/eJ* 


where  the  low  field  mobility,  fjQ,  exhibits  density  dependence  follow  the 
data  of  Blakemore  [14].  An  empirical  fit  to  this  data  yields 


Mo  =  5200  -  2800  tanh  (o . 902992 [ (logN)  -  16.845]} 


The  mobility  is  obtained  by  dividing  the  velocity,  Vn  by  the  electric  field, 
E.  For  holes, 


(i  ♦  (W) 


0W0 


where  |e|  is  the  magnitude  of  the  electric  field  and  Vsat  is  the 
saturated  drift  velocity. 

The  diffusivity  is  then  determined  using  the  Einstein  relation 


n  kT 
D  =  e~  M 
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where  k  is  Boltzmann's  constant,  and  T  is  the  temperature  (taken  here  as 
constant,  300°K) .  The  constants  used  in  the  mobility  expression  are  given 
in  Table  2. 


2.3  BOUNDARY  CONDITIONS  AND  DOPING  SPECIFICATION. 


The  specification  of  the  doping  distribution  and  the  boundary  conditions 
determine  the  type  of  device  under  consideration  and  the  bias  point.  For  the 
two-  and  three-dimensional  comparative  simulations  considered  here,  the  device 
structure  is  taken  as  two-dimensional  in  the  X-Z  plane.  Three-dimensional 
effects  are  introduced  by  the  presence  of  the  ionizing  particle  track.  Thus, a 
direct  evaluation  of  the  limitation  of  a  two-dimensional  modeling  approach  can 
be  made.  The  doping  distribution  for  the  device  is  thus  given  as 


C(x,z)  =  ND  -  Na 


Boundary  conditions  are  required  for  ohmic  contacts  and  free  surfaces.  The 
carrier  densities  at  ohmic  contacts  are  determined  through  the  assumption  of 
zero  space  charge 

N  -  P  =  C  (17) 


together  with  the  assumption  of  thermal  equilibrium 
NP  =  N£2 


The  concentrations  of  electrons  and  holes  at  ohmic  contacts,  as  well  as  the 
initial  distribution  of  carriers  throughout  the  device,  are  determined  by 
simultaneous  solution  of  Eqs .  (17  and  18). 

The  pof  •’ntial  at  ohmic  contacts  is  specified  relative  to  the  vacuum  level,  as 
it  must  be  if  a  consistent  treatment  of  heterojunctions  is  to  be  retained. 

At  N-type  contacts 


kT  N 
4  ”  VAPP  +  e~  ln  7T 


and  at  P-type  contacts 


kT  p  Ee  X 
V  ”  VApp  -  —  In  —  -  — ^  - 


It  should  be  noted  that  in  the  absence  of  any  variation  in  the  electron 
affinity,  and  under  the  assumption  that  the  Fermi  level  is  centered  between 


R 


e 


the  conduction  and  valence  bands  (i.e.  Nc  -  Nv) ,  Eq.  (19),  for  example, 
reduces  to  the  commonly  used  relationship  for  homojunctions , 


,  „  kT  ,  N 

*  "  vapp  +  T  ln  % 


(21) 


where  the  contributions  to  the  built-in  potential  from  the  band  gap  and 
electron  affinity  are  ignored  since  they  are  equal  at  all  contacts. 

At  free  surfaces,  the  normal  component  of  electric  field  and  current  density 
are  set  to  zero. 


2.4  MODELING  OF  THE  TRACK. 

Obviously,  modeling  of  the  track  in  three  dimensions  poses  no  significant 
problems.  The  generation  of  the  associated  excess  carriers,  on  a  per-unit 
volume  basis,  has  been  discussed  previously.  Thus,  all  that  need  be  specified 
in  three  dimensions  is  the  total  number  of  carriers  to  be  generated,  the  range 
of  the  ionizing  particle  (or  track  length)  and  an  assumption  of  the  initial 
diameter  of  the  track.  The  only  approximation  required,  when  using  a 
three-dimensional  Cartesian  coordinate  system,  is  that  the  initial  cross- 
section  of  the  track  be  modeled  by  a  series  of  square  cells.  Thus,  the 
initial  track  cross-section  will  only  approximate  a  circular  cross-section. 

The  accuracy  of  this  initial  approximation  will  depend  on  the  size  of  the  grid 
spacing  relative  to  the  track  diameter.  However,  even  when  the  grid  spacing 
is  on  the  order  of  the  track  radius,  this  initial  approximation  rapidly  decays 
to  a  cylindrical  distribution,  as  demonstrated  in  three-dimensional 
simulations  of  a  simple  diode  [8,9], 

In  two  dimensions,  however,  significant  approximations  must  be  made.  As 
discussed  in  [8],  in  two  dimensions  only  the  X-Z  plane  of  the  device  is 
considered,  and  the  excess  particle  generated  carriers  are  introduced  over  an 
area  of  this  plane.  The  length  of  this  region  is  taken  as  the  particle 
penetration  or  track  length  while  the  width  is  taken  as  the  initial  track 
diameter.  Both  the  device  and  the  track  extend  indefinitely  in  the  unmodeled 
third  dimension.  The  particle  track  is  thus  represented  by  a  slab  of  excess 
charge  rather  than  as  a  cylinder.  To  yield  meaningful  results  from  the 
simulations,  a  depth  in  the  third  dimension  must  be  specified.  When  choosing 
this  arbitrary  depth  several  constraints  provide  guidance.  However,  not  all 
of  the  constraints  can  be  met  simultaneously.  The  constraints  considered  are 
as  follows: 

(1)  the  total  ion -generated  excess  charge, 

(2)  the  density  of  carriers  in  the  track, 

(3)  the  volume  of  the  slab  representing  the  particle  track,  and 

(4)  the  physical  depth  of  the  device. 

Since  the  device  structure  is  itself  assumed  to  be  two-dimensional,  the  fourth 
constraint  also  sets  the  device  contact  areas. 
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It  is  obviously  necessary  to  introduce  the  correct  number  of  electron-hole 
pairs  for  a  given  particle,  thus  the  first  constraint  must  be  met.  The 
remaining  three  constraints  are  all  interdependent.  If  the  depth  of  the 
device  is  chosen  as  the  actual  physical  device  depth,  then  the  contact  areas 
will  be  correct  but  the  volume  of  the  slab  representing  the  particle  track 
will  be  much  greater  than  the  actual  cylindrical  track,  and  to  maintain  the 
proper  total  charge  the  track  density  must  be  reduced  below  its  actual  level. 
On  the  other  hand,  if  the  device  depth  is  chosen  to  be  on  the  order  of  the 
initial  track  diameter,  then  the  volume  of  the  slab  of  excess  carriers  will 
closely  approximate  the  actual  track  volume  and  the  carrier  density  within  the 
track  will  be  near  the  actual  value.  However,  the  device  contacts  will  be 
much  smaller  than  in  reality.  In  a  recent,  comparative  two-  and 
three-dimensional  study  performed  by  the  present  authors  [8] ,  it  was  found 
that  the  first  approach,  choosing  the  device  depth  as  the  relevant  parameter 
and  adjusting  the  track  density,  yielded  similar  current  pulses  for  both  two- 
and  three-dimensional  simulations;  the  second  choice  resulted  in  significantly 
larger  collection  times  and  over-predicted  field  distortions.  The  first 
approach  was  also  found  to  yield  excellent  results  in  the  NMOS  simulation 
reported  in  [10].  Thus,  the  approach  followed  in  the  two-dimensional  device 
simulation  performed  here  is  to  assign  a  value  to  the  depth  parameter 
representative  of  the  device  depth. 

It  should  be  cautioned,  however,  that  while  the  results  of  the  comparative 
studies  of  [8]  and  [10]  gave  good  predictions  of  the  current  pulse,  compared 
to  the  full  three-dimensional  simulation,  the  details  of  carrier  transport  in 
the  devices  were  significantly  different.  Therefore,  in  complex  devices,  such 
as  that  considered  here,  if  carrier  transport  details  are  important  the 
present  approach  may  be  only  qualitative  at  best.  For  this  reason,  the 
authors  continue  to  advocate  that  full  three-dimensional  simulations  be 
performed  for  complex  structures,  at  least  to  verify  the  two-dimensional 
modeling  approach.  This  comparison  and  verification  is  one  of  the  objectives 
of  the  present  research  applied  to  a  GaAs  JFET. 


SECTION  3 


THE  NUMERICAL  METHOD 


3.1  PHILOSOPHY  OF  THE  SOLUTION  PROCEDURE 

A  detailed  discussion  of  the  solution  technique,  including  the  development  of 
consistent  difference  approximations  to  the  governing  equations,  is  given  in 
[15],  Thus,  the  discussion  here  will  be  limited  to  the  philosophy  behind  the 
technique  and  implementation  on  a  vector  machine. 

In  an  effort  to  develop  a  highly  efficient  solution  technique  to  the  system  of 
Eqs.  (1,  2  and  11)  it  is  first  recognized  that  this  system  is  a  coupled, 
non-linear  system.  If  solved  as  a  coupled  system,  it  will  require  utilization 
of  methods  designed  for  coupled  elliptic  equations  to  obtain  a  solution  at 
each  time  step.  These  methods  typically  introduce  some  outer  iteration  to 
treat  nonlinearity  and  the  methods  used  to  solve  the  difference  approximations 
to  the  linearized  system  are  often  computationally  intensive  in  two 
dimensions,  and  totally  impractical  to  implement  in  three.  By  contrast,  the 
method  used  here  eliminates  nonlinear  iteration,  uses  a  noniterative  yet 
highly  efficient  procedure  to  solve  the  difference  approximations  to  the 
continuity  equations,  and  uses  an  extremely  efficient  iterative  technique  to 
solve  the  equation  governing  the  potential.  As  shall  be  discussed,  these 
solution  techniques  are  also  ideally  suited  for  implementation  on  vector 
machines  and/or  parallel  processors  making  them  even  more  attractive. 

The  first  step  in  this  procedure  is  to  decouple  Poisson's  equation  from  the 
continuity  equation  in  a  manner  which  does  not  adversely  effect  stability  of 
the  overall  solution  algorithm.  This  is  accomplished  by  reformulating  the 
continuity  equations  by  expanding  the  drift  term,  and  substituting  the  space 
charge  for  the  Laplacian  of  the  potential.  After  manipulation,  the  result  is 


“  -V  •  NMnVFn  -  VMn H  •  V0 


Ne 

Mn  l1  (N-P-C) 


+  7  Vf  •  +  V  •  DnVN  +  G  -  R 


al  -  v  •  P^p  VFp  +  VMpP  .  p  (N-P-C) 


-  r-p  I  Ve  •  Vtf.  +  V  •  DpVP  +  G  -  R 


To  ensure  conservation  of  total  current,  Poisson's  equation  is  recast  as  a 


(24) 


JIN. ,  1.1  M'tJ'l.l'l.t'IJ 


statement  of  total  current, 

—  *dt€V*  -  -  eV  •  N/in  +  Fn)  -eV  .  PMpV(tf>  +  Fp) 

+  eV  • (DnVN  -  DpVP) 


From  Eqs.  (22-24)  It  is  easily  shown  that  at  steady  state,  Poisson's  equation 
is  satisfied  exactly.  A  small  error  is  introduced  in  transients,  as  discussed 
in  [15], 

Eqs.  (22-24)  form  the  basis  of  the  solution  algorithm.  To  advance  the 
solution  from  tn  to  tn+^  -  tn  +  At,  the  mobilities  and  diffusivities 
are  evaluated  using  the  electric  field  and  carrier  densities  at  tn. 
Additionally,  the  gradients  of  Potential  appearing  in  the  second  and  fourth 
terms  on  the  R.H.S.  of  Eqs.  (22  and  23)  are  evaluated  at  the  tn  level.  This 
effectively  decouples  the  continuity  equations  from  the  total  current 
constraint  and  allows  the  carrier  concentrations  to  be  advanced  first.  This 
decoupling  does  not  introduce  a  stability  constraint  [15]. 

The  advance  in  time  of  the  carrier  concentrations  is  performed  by  solving  the 
continuity  equations  (Eqs.  (22  and  23))  as  a  block  2x2  coupled  system  through 
application  of  a  linearized  block  implicit  (LBI)  method  [16].  The  continuity 
equations  are  of  the  form 


D(*)  +  S(^) 


where  4>  =  (N,P)^,  D($)  represent  those  terms  in  Eqs.  (22  and  23)  which 
contain  spatial  derivatives  of  <f> ,  and  S(^)  represent  source  terms  such  as 
the  recombination,  generation  and  space  charge  terms.  Eq.  (25)  is  then  time 
differenced  using  a  backwards  differencing  scheme, 


D(^)n+1  +  S(^)n+1  +  0(At) 


where  A^n_r^-  =  -  <f>nt  and  the  superscripts  refer  to  the  time 

level.  D(<£)n+-*-  and  S(<f>)n+^~  are  then  formally  linearized  in  time  using 
a  Taylor  series  expansion  about  the  solution  at  time  tn  as 

G(^)n+1  =  G{4>)n  +  At  n_M_  n+  0(At2) 


Substituting  a  forward  difference  approximation  for  the  time  derivative  in 

Eq.  (26) 


G(^)n  +  Mil  %n+l  +  0(At2) 
o<p 
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Eq.  (26)  may  then  be  expressed  as 


(A  +  AtL)A^n+1  -  At 


^D(*)n 


+  S(*)n  +  O(At) 


where 


A  -  I  - 


At3^n 

d<t> 


t  _  3D (<f>)  n 

84> 


When  the  L  operator  is  approximated  by  three -point  difference  approximations 
Eq.  (29)  represents  a  block  2x2  matrix  equation  which  may  be  written  at  each 
grid  point  in  the  solution  domain.  In  one  dimension,  the  result  is  a  block 
2x2  tridiagonal  coefficient  matrix  which  may  be  solved  efficiently  using 
direct  block  tridiagonal  elimination.  For  two-  or  three-dimensional 
approximations,  while  the  block  size  remains  2x2,  since  it  is  determined  by 
the  number  of  coupled  equations ,  the  bandwidth  of  the  resulting  coefficient 
matrix  increases  significantly.  In  two  dimensions,  on  a  square  mesh  of  NxN 
points,  the  rank  of  the  coefficient  matrix  is  of  order  N2  and  the  bandwidth 
is  of  order  N.  In  three  dimensions,  with  an  NxNxN  mesh,  the  rank  of  the 
coefficient  matrix  is  of  order  N3  and  the  bandwidth  of  order  N2 . 

Obviously,  use  of  direct  inversion  techniques  for  such  matrices,  for  even 
relatively  small  meshes,  would  result  in  a  computationally  intensive  and 
inefficient  solution  procedure.  For  this  reason,  iterative  matrix  solvers  are 
often  used.  However,  for  large  meshes  even  these  iterative  solvers  may  become 
prohibitive  and  should  be  avoided,  if  possible.  This  can  be  accomplished,  due 
to  the  parabolic  nature  of  continuity  equations,  by  applying  a  consistently 
split  ADI  procedure  [17]  to  solve  Eq.  (29).  To  split  or  factor  Eq.  (29),  the 
L  operator  is  separated  into  its  directional  components,  L  =  1^  +  Ly  +  Lz, 
and  Eq.  (29)  is  rewritten  as  a  sequence  of  one -dimensional  systems  along  each 
mesh  line  in  the  x,  y  and  z  directions  respectively: 


(A  +  AtLx)A«^*  =  At|o(^)n  +  S  (<£)nj 


(A  +  AtLy)A^**  =  AAtf’1 


(A  +  AtLz)A^***  =  AA  4>** 


(32a) 


(32b) 


(32c) 


Elimination  of  the  intermediate  steps  in  Eq.  (32)  yields 

(A  +  AtLx)A’1(A  +  AtLy)A‘ 1  (A  +  AtLz)A^***  -  At|D(<£)n  +  S(<^)n] 


ft 
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Comparison  of  Eqs.  (29  and  33)  shows  that  Eq.  (33)  approximates  Eq.  (29)  to 
0(At2),  thus 

A^n+i  _  A^***  +  Q(At2)  (34) 


While  this  factorization  error  may  place  additional  restrictions  on  the  time 
step  when  considering  accuracy,  the  overall  reduction  in  computational  effort 
will  typically  more  than  offset  this  limitation.  Each  of  Eqs.  (32)  is 
tridiagonal,  thus  direct  elimination  can  be  implemented  in  a  highly  efficient 
manner.  On  an  NxNxN  mesh,  the  factorization  reduces  the  need  to  solve  a  rank 
N3  matrix  to  a  task  requiring  the  solution  of  3N2  tridiagonal  matrices  of 
rank  N.  Since  only  tridiagonal  matrices  need  to  be  solved,  regardless  of  the 
mesh  structure,  the  number  of  operations  per  mesh  point  remains  constant  and 
the  computational  effort  required  to  solve  the  continuity  equations  varies 
linearly  with  total  mesh  points  when  implemented  on  a  scalar  machine. 


Having  advanced  the  carrier  concentrations,  the  total  current  constraint,  Eq. 
(24)  must  be  solved  for  the  potential.  Since  the  carrier  densities  at  tn+1 
are  now  known,  Eq.  (24)  can  be  differenced  fully  implicitly  while  only 
remains  unknown.  It  must  be  noted  that  while  a  time  derivative  appears  in  Eq. 
(24) ,  this  equation  remains  elliptic  and  must  be  solved  iteratively.  To 
accomplish  this  Eq.  (24)  is  recast  as 


3V  • 


+  eV*Npn  V<V  +  Fn)  +  eV*P/ip V{ij>  +  Fp)  -  eV*(DnVN  -  DpVP) 


and  Eq.  (35)  expressed  as 


(A  +  iL)A^i+1  =  ijDW1  +  S(0)i] 


where 


T  _  i  dS< 

P  dxt> 


L  =  - 


mr> 


Here,  D(\6)  includes  the  first  three  terms  on  the  R.H.S.  of  Eq.  (35)  and 
S(V>)  the  last  term.  The  superscript  "i"  refers  to  an  iteration  index  and 
p  is  an  acceleration  parameter  which  varies  both  spatially  and  from 
iteration  to  iteration.  Eq .  (36)  may  be  ADI  split,  as  were  the  continuity 
equations,  following  Eqs.  (32).  However,  in  contrast  to  the  continuity 
equations,  the  total  current  constraint  must  be  iterated  to  convergence  at 
each  physical  time  step.  With  the  proper  choice  of  acceleration  parameters 
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this  can  be  accomplished  rapidly  and  efficiently.  It  should  also  be  noted 
that  since  the  D  operator  is  linear  in  L  -  -D  and  the  D  operator  need 
only  be  computed  at  the  start  of  the  iteration  and  stored.  Similarly,  the  S 
operator  is  not  a  function  of  ^  in  this  case,  thus  dS/difi  -  0.  The  S 
operator  also  need  be  computed  only  once  and  stored.  At  convergence, 

AV>i+1  -  A tf***  will  go  to  zero  and,  as  may  be  observed  from  Eq.  (33), 
the  factorization  error  will  also  go  to  zero  and  the  difference  approximations 
to  Eq.(24)  is  solved  exactly.  This  completes  the  advance  of  the  solution  from 
tn  to  tn+l.  The  process  is  then  repeated  for  the  next  time  step. 


3.2  IMPLEMENTATION  ON  A  VECTOR  MACHINE. 

The  procedure  described  in  the  previous  section  has  proven  to  be  a  very 
efficient  procedure.  The  basic  efficiency  of  the  algorithm  lies  in  the  low 
operations  count  of  the  procedures  used  to  solve  the  systems  of  difference 
equations.  On  a  three-dimensional  mesh  with  a  total  of  Nx  grid  points,  the 
number  of  operations  required  to  solve  Eqs.  (29  and  36)  is  given  as  [15] 


0PT3D  »  108NX  +  (15Nt) I 


where  I  is  the  number  of  iterations  required  in  the  solution  of  the  total 
current  constraint.  On  a  two-dimensional  mesh  this  reduces  to 


0PT2d  =  72NT  +  (10Nt)I 


since  only  two  of  the  ADI  sweeps  need  to  be  performed, 
problem  the  result  is 


For  a  one -dimensional 


0PT1D  *  36NT  +  (5Nt)I 


Note  that  in  each  case,  the  multi -dimensional  operation  count  is  simply  the 
one -dimensional  result  times  the  number  of  dimensions. 

In  three  dimensions,  the  total  number  of  grid  points  Nx  is  the  product  of 
the  number  of  mesh  points  in  each  direction,  respectively 


%  -  Nx  *  Ny  *  Nz 


Thus,  the  three-dimensional  operation  count  is  the  sum  of  the  operations 
required  for  each  of  Eqs.  (32a  -  32c) 


°PT3D  *  (36NX  +  5NXI)  *  NyNz 
+  (36Ny  +  5NyI)  *  NXNZ 

+  (36NZ  +  5NZI)  *  NxNy 
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where  there  are  NyNz  "x"  implicit  lines,  NXNZ  "y"  implicit  lines  and 
NxNy  "z"  implicit  lines.  Since  each  of  the  "x"  implicit  lines  requires 
the  same  sequence  of  operation  to  be  performed  to  solve  Eq.  (32a),  then  the 
solution  of  the  "x"  implicit  lines  can  be  done  in  parallel,  or  in  a  vector 
loop.  Similarly,  the  operations  required  to  set  up  or  generate  the  difference 
equations  which  constitute  Eq.  (32a)  can  also  be  performed  in  vector  loops. 

As  a  result,  rather  than  set  up  and  then  solve  Eq.  (32a)  along  each  of  the 
NyNz  "x"  implicit  lines  one  at  a  time,  all  NyNz  "x"  implicit  lines  may 
be  processed  at  once  in  a  single  vectored  loop,  providing  sufficient  storage 
is  available.  The  NyNz  scalar  operations  are  replaced  by  a  single  vector 
operation.  In  general,  a  vector  length,  VL,  may  be  defined,  and  the  number  of 
vector  loops  needed  to  solve  Eq.  (32a)  will  be  the  integral  part  of 
(NyNz)/VL.  Similar  results  are  obtained  for  the  "y"  and  "z"  implicit 
counterparts,  Eqs.  (32b  and  c).  The  number  of  vector  operations  required  to 
solve  Eq.  (32),  for  three-dimensional  problems  then  becomes 

OPTV3d  =  108NX/VL  +  (15NX/VL) I  (44) 


If  an  NxNxN  mesh  is  used  and  sufficient  storage  exists  for  VL  -  N2 
Nx2/3 ,  then 


OPTV3d  «  108Nt3  +  (15Nt3)I 


For  two-dimensional  problems  the  result  is 


0PTV2d  «  72Nt2  +  (10Nt2)I 


Comparison  of  Eqs.  (45  and  46)  with  their  scalar  counterparts,  Eqs.  (39  and 
40) ,  immediately  reveals  that  there  are  substantial  savings  for  both  two-  and 
three-dimensional  algorithms;  that  the  larger  the  problem,  the  greater  the 
potential  savings;  and  the  potential  savings  is  greater  in  the 
three-dimensional  case.  However,  in  the  present  simulations,  the  vector 
length  was  limited  to  VL  ■  Nx1/3  for  the  three-dimensional  problems 
and  VL  «  NXV2«N  for  two-dimensional  problems.  Thus,  while  Eq.  (46) 
yields  an  indication  of  the  performance  obtained  in  two  dimensions,  the 
three-dimensional  performance  followed  the  estimate 


0PTV3d  a  108Nt3  +  (15NX3)1 
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SECTION  4 


RESULTS 


4.1  DEVICE  STRUCTURE  CONSIDERED. 

The  structure  of  the  GaAs  JFET  device  considered  under  the  present  research  is 
shown  in  Fig.  1.  The  device  structure  is  taken  to  be  invariant  in  the  third 
dimension  (into  the  page)  thus  the  three-dimensional  effects  are  limited  to 
those  due  to  the  particle  track.  The  device,  which  is  representative  of 
current  JFET  technology,  consists  of  an  N^  source  and  drain,  doped  to 
2x10 17 /cm3,  a  p-gate  doped  to  2x10 17 /cm3  and  a  1.0x10 17 /cm3  N  type 
channel.  The  substrate  is  jt  GaAs,  nominally  p-type  (lxl014cm3).  The 
dimensions  of  the  device  are  shown  in  the  figure  but  it  is  worthy  to  note  that 
the  source-drain  spacing  is  3pm.  The  device  is  taken  as  25pm  wide  (into 
the  page).  The  source  and  drain  diffusions  penetrate  0.5pm  into  the 
device.  The  channel  is  0.3pm  deep  and  the  gate  diffusion  penetrates  0.2pm 
into  the  channel. 

The  device  was  biased  in  the  off  state  with  1  volt  applied  to  the  drain  and  0 
volts  applied  to  the  source  and  the  gate  contact.  The  substrate  was  also 
grounded . 

4.2  DETAILS  OF  THE  PARTICLE  TRACK. 

The  incident  particle  considered  in  the  present  simulations  is  a  5MeV  alpha 
particle.  As  shown  in  Fig.  1,  the  particle  is  assumed  to  enter  the  device 
through  the  gate  diffusion  and  penetrate  approximately  18.5pm,  the  range  of 
such  a  particle  in  GaAs  [18].  Since  4.5eV  are  required  to  produce  a  single 
electron-hole  pair  in  GaAs,  a  5MeV  particle  will  generate  about  9.375x10s 
electron-hole  pairs.  In  the  present  simulations  these  excess  carriers  are 
introduced  uniformly  along  the  track.  Since  the  initial  track  radius  is  of 
the  order  of  1000A,  this  results  in  a  track  density  of  1 . 38X101 8/cm3  for 
the  three-dimensional  simulation  where  the  track  is  represented  by  a  square 
cross-section  ,2x.2pm.  For  the  two-dimensional  simulation  the  track  is 
modeled  as  a  slab  ,2pm  thick,  18.5pm  long  and  extending  the  width  of  the 
device,  25pm.  Recall  that  in  two  dimensions  the  calculations  are  performed 
on  a  per-unit  width  basis  and,  as  discussed  in  section  2.4,  the  assignment  of 
the  assumed  width  of  the  device  is  critical.  Using  the  actual  device  width 
was  found  to  yield  similar  results  between  two-  and  three-dimensional 
simulations  in  [8]  and  [10]  and  this  approach  is  followed  here.  Using  25pm 
as  the  width  parameter  yields  a  track  density  of  1 . 103X101 6/cm3  substan¬ 
tially  below  the  actual  track  density,  but  yielding  the  correct  total  number 
of  excess  carriers. 


4.3  THE  COMPUTATIONAL  PROCEDURE. 

The  first  step  in  performing  the  simulations  is  the  construction  of  a  suitable 
finite  difference  grid.  For  the  two-dimensional  simulation  an  unequally- 
spaced  mesh  of  51  x  points  by  44  z  points  was  used.  The  mesh  in  the  X-Z  plane 


is  shown  in  Fig.  2.  The  mesh  is  structured  to  give  adequate  resolution  of 
junctions  and  the  particle  track.  For  the  three-dimensional  simulation,  the 
same  mesh  was  used  in  the  X-Z  plane  with  18  grid  points .unequally  spaced,  in 
the  y  direction.  The  X-Y  and  Y-Z  mesh  structures  are  shown  in  Fig.  3.  Note 
that  since  the  particle  is  assumed  to  strike  the  center  of  the  device,  the 
solution  will  be  symmetric  about  the  center,  and  only  half  the  device  need  be 
modeled  in  three  dimensions.  Thus,  the  y  grid  extends  from  the  center  of  the 
device  (and  particle  track)  to  the  device  edge.  Note  the  orientation  of  the 
coordinates  in  Figs.  1  through  3  as  this  will  aid  in  the  interpretation  of 
contour  plots  presented  later. 

The  next  step  is  to  compute  the  solution  for  the  initial,  undisturbed  state. 
Since,  in  the  absence  of  the  particle  track,  this  solution  is  invariant  in  the 
y  direction,  the  solution  need  only  be  computed  in  two  dimensions.  This 
result  was  used  for  the  initial  condition  in  both  the  two-dimensional  and 
three-dimensional  transient  charge  collection  simulations. 

The  two-dimensional  transient  simulation  was  then  computed  using  2000  time 
steps  and  a  CPU  time  on  a  Cray  X/MP24  computer  of  14.13  minutes  or  0.00018 
sec/grid  point/time  step.  This  was  14  times  faster  than  the  original  scalar 
version  of  the  code.  In  the  three-dimensional  simulation  smaller  time  steps 
were  required  initially  to  maintain  sufficient  accuracy  with  the  higher 
initial  track  density.  Thus,  the  three-dimensional  simulation  required  2075 
time  steps  yielding  a  CPU  time,  on  a  Cray  X/MP24,  of  6.83  hours  or  about 
0.00029  sec/grid  point/time  step.  This  was  roughly  15  times  faster  than  the 
projected  run  time  based  on  the  original,  scalar  code. 

4.4  THE  EQUILIBRIUM  SOLUTION. 

Since  the  equilibrium  solution  is  not  a  function  of  the  Y  coordinate 
direction,  it  need  be  examined  in  only  a  single  X-Z  plane.  In  Fig.  4  contours 
of  the  log  of  electron  and  hole  concentrations,  and  potential  are  shown.  The 
density  contours  are  equally  spaced  with  an  increment  one  order  of  magnitude. 
The  potential  contours  are  also  equally  spaced  with  an  increment  of  0.1 
volts.  The  results  are  as  expected,  showing  the  N+  diffusions  and  depletion 
of  the  electrons  in  the  channel.  Figs.  5  and  6  show  details  of  the  contours 
of  the  carrier  densities  in  the  channel  region  of  the  device  in  both  a  linear 
(Figs.  5a  and  6a)  and  a  log  scale.  The  tick  mark  on  the  bottom  of  the  figures 
is  2.5pm  from  the  contact  surface.  Fig.  7  shows  an  enlargement  of  the 
potential  contours  in  the  same  region  and  a  surface  plot  of  the  steady  state 
potential  distribution.  The  typical  high  field  region  at  the  drain  side  of 
the  gate  is  clearly  evident  in  this  figure. 

The  device  structure  is  such  that  it  is  a  normally  off  device,  thus  at  the 
applied  bias  level,  there  is  negligible  source-drain  current.  This  provides 
an  excellent  initial  condition  for  the  charge  collection  simulations  to  be 
performed.  Any  current  observed  during  the  transient  will  be  a  direct  result 
of  the  disturbance  generated  by  the  particle  strike. 

4.5  COMPARISON  OF  THE  TWO-  AND  THREE-DIMENSIONAL  TRANSIENT  RESULTS 

The  main  point  of  interest  in  the  comparison  of  the  two-  and  three-dimensional 


charge  collection  transient  simulations  is  the  shape  and  duration  of  the 
current  pulse  at  the  struck  gate  node  and  the  charge  collected  there.  The 
current  pulses  at  other  contacts  are  also  of  interest  as  they  may  affect 
device  operation.  The  current  pulses  at  the  source,  gate  and  drain  contacts 
are  shown  in  Fig.  8  for  the  two-dimensional  case.  The  positive  drain  current 
indicates  electrons  leaving  the  device  while  the  negative  source  current 
represents  injection  of  electrons.  The  gate  current  is  made  up  of  electron, 
hole  and  displacement  current,  however,  it  is  predominantly  hole  current  and 
represents  collection  of  holes  from  the  particle  track.  The  results  shown  in 
Fig.  8  are  very  similar,  both  in  duration  and  magnitude,  to  the  results 
obtained  for  a  JFET  in  [4]  despite  the  fact  that  the  source,  drain  and  gate 
diffusions  were  more  heavily  doped  (5xl018/cm3)  in  [4]. 

As  shown  in  Fig.  8,  the  gate  current  in  the  two-dimensional  simulation 
increases  rapidly  to  a  peak  of  0.325ma  in  less  than  lOpsec  and  then  decays  in 
an  exponential  manner.  By  lOOpsec  the  gate  current  is  approximately  0.02ma 
and  further  decay  takes  place  on  a  very  long  (relatively)  time  scale.  By 
400psec  the  gate  current  is  only  about  0.006ma  and  continues  to  decay  slowly. 

The  initial  peak  in  the  gate  current  is  associated  with  the  collection  of 
holes  generated  in  the  gate  depletion  region.  Once  these  holes  are  collected 
the  gate  current  decays  and  additional  holes  are  collected  slowly  as  they 
diffuse  into  the  junction  region.  As  will  be  discussed  subsequently,  field 
funneling  at  the  P+N  junction  between  the  gate  and  channel  does  not  occur, 
at  least  not  in  the  usual  sense. 

The  current  pulses  at  the  source  and  drain  are  also  shown  in  Fig.  8.  The 
drain  current  exhibits  a  rapid  rise  to  its  initial  peak  which  occurs  at 
approximately  lOpsec.  This  initial  peak  is  associated  with  the  collection,  of 
the  electrons  generated  in  the  gate  junction.  A  small  positive  peak  is  also 
observed  at  the  source.  This  is  also  associated  with  the  collection  of 
electrons  generated  in  the  gate  junction.  However,  this  peak  rapidly  decays 
and  the  source  current  becomes  negative,  balancing  the  drain  current,  which  is 
well  in  excess  of  the  gate  current.  Obviously  the  source-drain  current  is  not 
simply  associated  with  the  collection  of  the  excess,  particle  generated 
carriers.  Rather,  it  is  a  result  of  the  spreading  of  the  excess  carriers 
generated  in  the  substrate.  As  in  the  simulation  reported  in  [4],  the  excess 
carriers  lower  the  substrate  resistance  and  open  a  current  path  around  the 
gate  allowing  a  source-drain  current  to  flow.  This  current  will  continue  to 
flow  for  a  substantially  longer  time  than  that  associated  with  the  gate 
transient  since  the  carriers  in  the  substrate  move  primarily  by  diffusion  and 
are  collected  slowly,  or  recombined.  Since,  as  time  goes  on,  the  excess 
carrier  concentrations  become  small,  the  diffusion  and  recombination  processes 
will  become  even  slower  and  the  elimination  of  excess  carriers  from  the 
substrate  will  drag  on.  The  source -drain  current  will  continued  to  flow 
during  much  of  this  time. 

The  current  pulses  at  the  source,  gate  and  drain  contacts  for  the 
three-dimensional  simulation  are  shown  in  Fig.  9.  Here  it  is  observed  that 
while  the  results  are  qualitatively  similar,  quantitative  differences  are 
apparent.  The  peak  in  the  gate  current  is  only  0.18ma.  However,  after  25psec 
the  gate  currents  are  almost  identical.  The  drain  current  does  not  exhibit 
the  rapid  rise  to  an  initial  peak  and  the  positive  peak  in  the  source  current 


is  greater.  Additionally,  the  source-drain  current  is  not  as  high. 


The  charge  collected  at  the  gate  contact  for  the  two-  and  three-dimensional 
cases  is  shown  in  Figs.  10  and  11,  respectively.  As  shown,  the  results 
reflect  the  differences  in  the  initial  peak  in  the  gate  current.  After 
400psec  the  two-dimensional  result  indicates  the  collection  of  more  than  6%  of 
the  particle-generated  charge,  or  charge  contained  in  the  first  1.1pm  of  the 
track.  By  comparison  the  three-dimensional  result  shows  only  slightly  more 
than  5%  of  the  charge  collected,  corresponding  to  the  first  0.925pm  of  the 
track.  By  400psec,  however,  the  rate  of  charge  collection  is  nearly  the  same 
for  both  cases. 

The  differences  in  the  results  for  the  two-  and  three-dimensional  cases  are 
somewhat  surprising  in  light  of  the  excellent  agreement  obtained  in  [8]  for  a 
silicon  diode,  and  [10]  for  an  NMOS  device  using  a  similar  method  to  scale  the 
particle  track.  However,  the  discrepancy  between  the  present  results  and 
those  of  [8]  and  [10]  can  be  explained  by  consideration  of  differences  in  the 
device  structures  and  the  response  of  the  field  to  the  particle  track.  In  [8] 
and  [10]  the  diffusions  of  the  struck  N*  regions  were  deep  enough  so  that  at 
the  contacts  the  electric  field  was  very  small.  Since  the  carrier  density  at 
the  contact  was  fixed,  and  since  the  low  field  mobility  in  the  region  of  the 
contact  was  also  constant,  the  contact  resistance  was  also  constant  giving 
rise  to  a  linear  current-field  relationship  at  the  contact.  Thus,  in  the 
scaled  two-dimensional  simulations  of  [8]  and  [10],  while  the  field 
disturbance  at  the  contact  was  smaller  than  the  localized  disturbance  in  the 
contact  field  observed  in  the  three-dimensional  results,  the  two-dimensional 
disturbance  occurred  over  the  entire  device  width.  As  discussed  in  [8]  the 
result  is  that  the  average  field  over  the  contact  area  is  the  same  in  both 
cases  and,  since  a  linear  current-field  relationship  exists,  the  contact 
currents  and  collected  charge  were  equal.  In  the  present  JFET  simulations  the 
gate  P+  diffusion  is  only  on  the  order  of  2000A  deep  and  space  charge 
effects  are  present  near  the  contact  surface.  As  a  result  the  field  at  the 
contact  is  not  low.  In  fact  it  is  high  enough  to  result  in  velocity 
saturation.  Therefore,  a  linear  current-field  relationship  does  not  exist  and 
even  with  an  equal  average  field  over  the  contact  area  the  contact  currents 
are  not  equal.  The  higher  local  field  in  the  three-dimensional  case  does  not 
give  rise  to  a  linear  increase  in  the  local  current  density  because  the 
velocity  is  in  saturation.  Rather,  the  increase  in  the  field  results  in  a 
lower  diffusivity  which  upsets  the  local  balance  between  drift  and  diffusion 
and  modulates  the  local  current  density  in  a  nonlinear  fashion.  In  the 
present  case  this  results  in  a  lower  initial  gate  current  for  the 
three-dimensional  simulation. 

The  localization  of  the  field  distortions  in  the  three-dimensional  case,  and 
the  absence  of  the  field  tunneling  effect  at  the  gate  diffusion  are 
illustrated  in  Fig.  12.  Here  the  potential  along  the  axis  of  the  particle 
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track  is  compared  for  the  two-  and  three-dimensional  cases  with  the  steady 
state  result.  The  potential  along  a  line  parallel  to  the  track  axis  in  an  X-Z 
plane  12.5pm  below  the  track,  on  the  surface  of  the  device,  is  also  shown. 

The  results  are  at  40psec  following  the  particle  strike.  The  gate  contact  is 
at  a  distance  of  zero.  The  potential  distribution  along  the  track  axis  shows 
reasonable  agreement  between  the  scaled  two-dimensional  simulation  and  the 
three-dimensional  result.  It  is  readily  apparent  that  there  is  no  field 
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funneling  effect  down  the  axis  of  the  track,  drawing  holes  toward  the  gate. 
Instead,  the  field  in  this  region  increases  significantly,  with  the  greater 
increase  occurring  in  the  three-dimensional  case.  However,  the  results  can  be 
interpreted  as  a  field  funneling  effect  which  draws  electrons  from  both  the 
gate  diffusion  and  the  ir-type  substrate  into  the  N  channel.  The  comparison 
12.5/im  below  the  track  axis  shows  that  in  the  gate-channel  junction  there  is 
very  little  effect  of  the  disturbance  although  the  disturbance  is  felt 
slightly  in  the  substrate. 


The  distortion  of  the  potential  along  the  axis  of  the  track  and  elsewhere 
within  the  device  can  be  explained  by  examining  the  net  electron  and  hole 
charge  collected  at  all  device  contacts.  This  result  for  the  two-dimensional 
simulation  is  shown  in  Fig.  13  and  for  the  three-dimensional  simulation  in 
Fig.  14.  As  can  be  seen,  during  the  initial  phase  of  the  transient  more 
electrons  are  collected  than  holes.  Thus,  the  net  excess  charge  in  the  device 
is  positive,  and  in  accordance  with  Poisson's  equation  and  Green's  theorem  the 
curvature  of  the  potential  surface  must  become  more  negative,  and  the  field  at 
the  device  contacts  must  change  in  a  consistent  manner.  This  is  discussed  in 
detail  for  the  case  of  a  simple  two -terminal  diode  in  [11]  and  the  arguments 
brought  forth  there  apply  equally  to  any  device.  The  distorted  potential 
along  the  track  axis  is  consistent  with  these  requirements. 


It  should  be  observed  in  Figs.  13  and  14  that  after  approximately  150psec,  the 
net  charge  imbalance  remains  constant.  This  imbalance  is  expected  to  remain 
nearly  constant  until  all  the  excess  electrons  are  collected  from  the  device 
or  recombined.  At  that  time  the  source-drain  current  will  decay  to  zero  and 
the  remaining  excess  holes  will  be  collected  (at  the  gate  and  substrate)  and 
the  potential  will  return  to  its  original  state.  Based  on  a  constant 
collection  rate  for  holes  obtained  from  the  slope  of  the  curve  in  Fig.  14  at 
400psec,  the  time  elapsed  before  complete  field  restoration  (complete  charge 
collection)  will  be  at  least  another  17.4nsec.  Since  the  gate  current,  where 
most  of  the  holes  are  collected,  continues  to  decrease,  as  shown  in  Figs.  8 
and  9,  the  average  collection  rate  will  typically  be  at  most  half  the  value  at 
400psec ,  and  the  time  for  complete  restoration  may  be  in  excess  of  35nsec. 

This  time  scale  may  be  reduced  by  recombination,  however,  as  the  carriers 
diffuse  throughout  the  substrate  the  excess  carrier  densities  become  very  low 
and  the  recombination  effects  will  become  unimportant. 


Fig.  15  shows  contours  of  the  log  of  the  electron  concentration  at  five 
instants  during  the  transient  for  the  two-dimensional  case.  In  these  figures, 
as  in  Fig.  4,  the  contour  spacing  is  one  order  of  magnitude.  The  contacts  are 
on  the  left  side  of  the  figure  with  the  source  contact  at  the  bottom  of  the 
left  side.  In  Figs.  15a,  b  and  c  the  contour  value  in  the  center  of  the  track 
is  2xl0ls/cm3.  In  Figs.  15d  and  e,  only  the  2xl014/cm3  contour 
remains.  Similar  results  are  shown  in  Fig.  16  for  the  X-Z  plane  containing 
the  particle  track  in  the  three-dimensional  simulation.  Again,  the  contour 
spacing  is  one  order  of  magnitude,  however,  the  maximum  contour  values  at  the 
center  of  the  track  are  larger.  In  Fig.  16a  the  center  contour  is  at 
2xl017/cm3.  The  center  contour  value  in  Figs.  16b,  c  and  d  is 
2x10* 6/cm3  and  in  Fig.  16d  the  value  is  2xl016/cm3.  Thus,  while  the 
spreading  of  the  excess  charge  appears  qualitatively  similar  in  the  two-  and 
three-dimensional  simulations  quantitatively  they  are  quite  different.  This 
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Is  a  direct  result  of  the  higher  initial  track  density  of  the 
three-dimensional  simulation  and  shows  the  stronger,  more  localized 
disturbance  present  in  three  dimensions.  Also  note  the  difference  in  the 
contours  in  the  immediate  region  of  the  gate  diffusion. 

Similar  contours  for  hole  concentrations  are  shown  in  Figs.  17  and  18.  In 
Fig.  17,  for  the  two-dimensional  results,  the  maximum  contour  value  associated 
with  the  track  is  2xlOI6/cm3  in  17a  through  d  and  2xl014/cms  in  17e. 

In  Fig.  18  the  maximum  contour  values  are  2xlQ17/cm3  for  18a, 

2xl016/cm3  for  18b,  c  and  d,  and  2xl015/cm3  for  18e. 

Potential  contours  for  the  two-  and  three-dimensional  simulations,  in  the  X-Z 
plane  of  the  particle  track  are  shown  in  Figs.  19  and  20  at  the  same  instants 
in  time  as  the  electron  and  hole  contours.  In  these  figures  the  increment 
between  potential  contours  is  0.1  volts.  There  is  good  quantitative  and 
qualitative  agreement  between  these  comparisons  indicating  that  the  scaling  of 
the  particle  track  density  still  holds  even  though  the  limited  gate  current 
alters  the  current  and  charge  collection  transient  somewhat.  The  slightly 
greater  initial  disturbance  of  the  potential  along  the  track  axis  is  clearly 
evident  in  Fig.  20.  A  more  vivid  comparison  of  the  potential  disturbance  in 
this  plane  is  shown  in  Figs.  21  through  25  where  the  2-D  and  3-D  potential 
surfaces  are  compared  side-by-side.  The  similarity  in  the  distortion  of  the 
potential  surfaces  in  the  substrate  is  clearly  evident  here. 

While  the  comparisons  of  the  potential  surface  in  the  plane  of  the  particle 
track  presented  in  Figs.  21  through  25  show  both  qualitative  and  quantitative 
agreement  it  must  be  recalled  that  the  distortions  associated  with  the 
two-dimensional  result  extend  over  the  entire  assumed  device  width.  However, 
the  three-dimensional  disturbance  may  exhibit  significant  variation  over  the 
device  width.  To  examine  how  far  the  disturbance  caused  by  particle  track 
extends  in  the  third  (y)  dimension,  the  potential  variation  in  two  planes 
normal  to  the  particle  track  is  presented  in  Figs.  26  through  30  at  five 
instants  during  the  transient.  In  these  figures  the  viewer  is  looking 
directly  down  the  particle  track,  through  the  gate  contact.  Since  the 
solution  is  symmetric  about  the  z  axis,  only  the  top  half  of  the  device  is 
shown  in  these  figures.  Thus,  the  particle  track  is  centered  at  the  bottom  of 
the  figures.  In  each  figure,  the  left-hand  frame  represents  the  Y-Z  plane 
0.5/im  below  the  contact  surface,  just  at  the  edge  of  the  channel.  The 
right-hand  frame  is  for  a  Y-Z  plane  8.24/jm  below  the  contact  surface,  in  the 
substrate.  The  source  contact  is  on  the  right,  and  the  drain  is  at  the  left 
of  each  frame.  Note  that  in  the  channel  region  of  the  device  the  cisturbance 
never  extends  much  further  than  2 . 5^m  outward  from  the  track  axis  in  the  y 
direction.  The  nearly  parallel,  vertical  potential  contours  indicate 
two-dimensionality,  i.e.,  a  Tack  of  variation  in  the  y  direction.  In  the 
substrate,  however,  the  result  is  somewhat  different.  In  the  absence  of  the 
particle  track  this  plane  would  be  at  a  constant  potential.  The  particle  track 
creates  a  disturbance  which  propagates  radially  outward  from  the  track  axis. 
Since  there  is  no  spatial  variation  of  the  donor  or  acceptor  concentrations  in 
the  substrate,  the  disturbance  is  cylindrical  about  the  track  axis  as  opposed 
to  the  three-dimensional  behavior  observed  in  the  active  region  of  the  device. 
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SECTION  5 


SUMMARY  AND  CONCLUSIONS 


The  results  of  the  present  study  show  that  the  response  of  a  typical  GaAs  JFET 
structure  to  single  particle  radiation  entails  a  much  more  complex  process 
than  simple  field  funneling  at  the  struck  node.  In  the  present  study  a  JFET 
device,  which  is  normally  in  the  off  state,  was  struck  through  the  gate 
diffusion  at  normal  incidence.  However,  the  collection  of  charge  at  the  gate 
was  not  augmented  by  the  field  funneling  effect.  Instead,  a  field  funneling 
effect  was  set  up  between  the  gate  and  the  N-channel,  and  the  channel  and  the 
substrate  drawing  excess  electrons  into  the  channel  to  subsequently  be 
collected  at  the  drain.  This  was  generated  by  a  charge  imbalance  which 
resulted  from  the  rapid  collection  of  electrons  from  the  gate  depletion  region 
while  the  holes,  due  to  their  lower  transport  coefficients  and  velocity 
saturation  at  the  gate  contact,  were  collected  more  slowly.  The  charge 
imbalance  resulted  in  the  distortion  of  the  potential  surface  in  an  attempt  to 
maintain  charge  neutrality.  It  appeared  that  this  initial  imbalance  would 
remain  constant  for  the  duration  of  the  event.  However,  a  more  significant 
aspect  of  the  response  was  that  the  spreading  of  the  excess  charge  throughout 
the  substrate  lowered  the  substrate  resistance  significantly  and  opened  a 
source -drain  current  path  in  the  off  device.  In  essence,  the  presence  of  the 
N-channel  blocked  the  collection  of  holes  from  the  substrate  and,  in  order  to 
maintain  a  nearly  charge-neutral  state,  electrons  were  injected  at  the  source 
to  replace  those  collected  at  the  drain.  There  were,  in  effect,  two 
components  associated  with  the  drain  current.  First,  there  was  a  small 
component  of  the  drain  current,  balanced  by  a  hole  current  at  the  gate,  which 
accounted  for  collection  of  the  excess  carriers .  As  the  results  indicated,  by 
AOOpsec  this  fraction  of  the  drain  current  was  very  small.  Second,  there  was 
the  much  larger  portion  of  the  drain  current  balanced  by  the  source.  This  was 
similar  to  a  typical  FET  current  and  was  due  to  a  modulation  of  the  field  in 
the  channel  and  substrate  by  the  spreading  of  the  excess  carriers.  The 
simulations  indicated  that  this  would  continue  until  the  excess  holes  were 
collected.  Since  the  excess  carriers  move  through  the  substrate  primarily  by 
diffusion  and  since  the  holes  must  diffuse  across  the  channel  to  the  gate  to 
be  collected,  the  source-drain  current  flowed  for  a  much  longer  time  than  that 
associated  with  the  gate  transient.  The  effect  is  exaggerated  in  GaAs  due  to 
the  low  diffusivity  of  holes.  The  present  simulations  indicated  that  this 
source-drain  current  could  continue  to  flow  35nsec  after  the  particle  strike 
even  though  the  gate  transient  was  effectively  over  in  400psec . 


In  two  previous  studies  [8,10]  comparative  two-  and  three-dimensional 
simulations  of  a  silicon  diode  and  an  NMOS  structure  demonstrated  that 
appropriate  scaling  of  the  charge  density  along  an  ion  track  in  a 
two-dimensional  simulation  yielded  a  current  pulse  and  charge  collection  rate 
that  was  in  excellent  qualitative  and  quantitative  agreement  with  the 
three-dimensional  results.  Thus,  it  appeared  that  SEU  simulations  of  such 
devices  could  be  performed  with  reasonable  confidence  using  the  two-dimen¬ 
sional  approach.  The  present  study  was  undertaken  to  further  verify  such  a 
scaling  approach  in  a  GaAs  JFET.  While  the  results  of  the  present  simulations 
showed  qualitative  agreement  between  the  two-  and  three-dimensional  JFET 
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simulations,  quantitative  differences  in  the  current  pulse  at  the  device 
contacts  and  the  resulting  charge  collection  rate  were  observed.  The  peak 
current  at  the  struck  gate  node  was  significantly  lower  in  the 
three-dimensional  simulation,  by  almost  a  factor  of  two.  After  25psec 
however,  the  gate  currents  were  almost  identical  in  both  cases.  After 
400psec,  the  gate  current  had  decayed  to  less  than  4%  of  the  peak  value.  At 
that  time  the  collected  charge  at  the  gate  node  was  about  6%  of  the  total  ion 
generated  charge  in  the  two-dimensional  case  and  about  5%  in  the 
three-dimensional  case  due  to  the  lower  peak  current. 

Significant  differences  were  also  observed  in  the  resulting  source-drain 
current  which  flows  due  to  a  lowering  of  the  substrate  resistance  by  the 
spreading  of  the  excess  carriers  generated  there.  Here,  the  drain  current  did 
not  rise  to  as  high  an  initial  peak  in  the  three-dimensional  simulation  and 
the  long  term  source -drain  current  level  was  somewhat  lower. 

After  careful  examination  and  comparison  of  the  simulation  results  it  was 
determined  that  the  primary  cause  of  the  differences  in  the  two-  and 
three-dimensional  results  was  a  result  of  the  device  structure.  Because  of 
the  shallowness  of  the  p+  gate  diffusion,  velocity  saturation  occurred  at 
this  contact,  even  before  the  particle  strike,  and  the  contact  resistance  is 
not  constant  but  highly  dependent,  in  a  nonlinear  fashion,  on  the  field  at 
this  contact.  In  the  two-dimensional  case,  the  field  distortions  at  this 
contact  are  uniform  across  the  assumed  device  depth.  However,  in  the 
three-dimensional  case,  only  localized  field  distortions  occur  and  the  average 
contact  resistance  is  different  from  the  two-dimensional  case,  resulting  in 
the  observed  difference  in  the  current  pulse  at  the  gate  contact.  In  the 
results  reported  in  [8,10]  the  contact  resistances  were  constant  at  all 
contacts.  Thus  it  is  concluded  here  that  scaling  of  the  track  density  in 
two-dimensional  simulations  will  provide  accurate  results  for  current  pulses 
and  charge  collection  transients  only  if  the  contact  resistances  are 
constant.  This,  in  effect,  requires  that  the  field  at  the  contacts  remain  low 
enough  to  maintain  operation  in  the  constant  mobility  region  of  the  velocity 
field  relationship.  Outside  this  range,  while  qualitatively  accurate  results 
may  be  obtained,  the  magnitude  of  any  quantitative  differences  can  only  be 
accessed  through  accurate,  three-dimensional  simulations.  It  is  therefore 
apparent  that  in  the  design  of  devices  which  are  hardened  against  single 
particle  radiation  effects,  two-dimensional  simulation  can  provide  a  rapid  and 
efficient  means  of  exploring  various  device  modif ications  and  structural 
variations.  However,  once  such  a  screening  process  is  completed,  the  results 
of  effective  hardening  approaches  should  be  verified,  at  least  selectively, 
using  three-dimensional  simulation  where  possible. 
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TABLE  1 


AVALANCHE  GENERATION  COEFFICIENTS 


1.90  xIO5 
2.21  xIO5 
5.75  x10s 
6.57  x10s 
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Figure  6.  Steady  state  contours  of  hole  density  in  the  channel 

region  of  the  JFET.  a)  contours  of  P,  b)  contours  of  log  P 


Current  pulses  as  a  function  of  time  at  source,  drain 
and  gate  contacts  after  particle  strike,  2-D  result. 
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Figure  12.  Comparison  of  potential  distribution  along  the  track  axis 
for  the  2-D  and  3-D  simulations  with  the  steady  state 
result  and  the  3-D  result  12.5pm  blow  the  track  axis 
40  psec  following  the  particle  strike. 


Figure  15.  Contoura  of  log  N  at  a)  10  paec,  b)  40  paac,  c)  100  paec, 
d)  200  paec,  and  a)  400  pa«c  following  tha  particle  strike 
for  tha  2-D  case. 


•WV 


Distribution  list  for  "Simulation  of  the  Response  of  a  Gallium  Arsenide  JFET 
to  Single  Particle  Radiation  in  Two  and  Three  Dimensions" . 


No.  Copies 


Naval  Research  Laboratory  1 

Attn:  Dr.  Edward  Peterson,  Code  4611 
4555  Overlook  Avenue,  S.W. 

Washington,  DC  20375 


Mr.  LeRoy  Dombkowski  1 

UAA/J9 

DCASMA 

96  Murphy  Road 
Hartford,  CT  06114-2173 


Director  6 

Naval  Research  Laboratory 
Attn:  Code  2627 
Washington,  DC  20375 


Defense  Technology  Information  Center 
Bldg.  5,  Cameron  Station 
Alexandria,  VA  22314 


12 


